home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_k.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.4 KB  |  244 lines
  1. /* specfunc/bessel_k.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_pow_int.h>
  26. #include <gsl/gsl_sf_gamma.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28.  
  29. #include "error.h"
  30. #include "check.h"
  31.  
  32. #include "bessel.h"
  33.  
  34. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  35.  
  36. /* [Abramowitz+Stegun, 10.2.4 + 10.2.6]
  37.  * with lmax=15, precision ~ 15D for x < 3
  38.  *
  39.  * assumes l >= 1
  40.  */
  41. static int bessel_kl_scaled_small_x(int l, const double x, gsl_sf_result * result)
  42. {
  43.   gsl_sf_result num_fact;
  44.   double den  = gsl_sf_pow_int(x, l+1);
  45.   int stat_df = gsl_sf_doublefact_e((unsigned int) (2*l-1), &num_fact);
  46.  
  47.   if(stat_df != GSL_SUCCESS || den == 0.0) {
  48.     OVERFLOW_ERROR(result);
  49.   }
  50.   else {
  51.     const int lmax = 50;
  52.     gsl_sf_result ipos_term;
  53.     double ineg_term;
  54.     double sgn = (GSL_IS_ODD(l) ? -1.0 : 1.0);
  55.     double ex  = exp(x);
  56.     double t = 0.5*x*x;
  57.     double sum = 1.0;
  58.     double t_coeff = 1.0;
  59.     double t_power = 1.0;
  60.     double delta;
  61.     int stat_il;
  62.     int i;
  63.  
  64.     for(i=1; i<lmax; i++) {
  65.       t_coeff /= i*(2*(i-l) - 1);
  66.       t_power *= t;
  67.       delta = t_power*t_coeff;
  68.       sum += delta;
  69.       if(fabs(delta/sum) < GSL_DBL_EPSILON) break;
  70.     }
  71.  
  72.     stat_il = gsl_sf_bessel_il_scaled_e(l, x, &ipos_term);
  73.     ineg_term =  sgn * num_fact.val/den * sum;
  74.     result->val = -sgn * 0.5*M_PI * (ex*ipos_term.val - ineg_term);
  75.     result->val *= ex;
  76.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  77.     return stat_il;
  78.   }
  79. }
  80.  
  81.  
  82. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  83.  
  84. int gsl_sf_bessel_k0_scaled_e(const double x, gsl_sf_result * result)
  85. {
  86.   /* CHECK_POINTER(result) */
  87.  
  88.   if(x <= 0.0) {
  89.     DOMAIN_ERROR(result);
  90.   }
  91.   else {
  92.     result->val = M_PI/(2.0*x);
  93.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  94.     CHECK_UNDERFLOW(result);
  95.     return GSL_SUCCESS;
  96.   }
  97. }
  98.  
  99.  
  100. int gsl_sf_bessel_k1_scaled_e(const double x, gsl_sf_result * result)
  101. {
  102.   /* CHECK_POINTER(result) */
  103.  
  104.   if(x <= 0.0) {
  105.     DOMAIN_ERROR(result);
  106.   }
  107.   else if(x < (M_SQRTPI+1.0)/(M_SQRT2*GSL_SQRT_DBL_MAX)) {
  108.     OVERFLOW_ERROR(result);
  109.   }
  110.   else {
  111.     result->val = M_PI/(2.0*x) * (1.0 + 1.0/x);
  112.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  113.     CHECK_UNDERFLOW(result);
  114.     return GSL_SUCCESS;
  115.   }
  116. }
  117.  
  118.  
  119. int gsl_sf_bessel_k2_scaled_e(const double x, gsl_sf_result * result)
  120. {
  121.   /* CHECK_POINTER(result) */
  122.  
  123.   if(x <= 0.0) {
  124.     DOMAIN_ERROR(result);
  125.   }
  126.   else if(x < 2.0/GSL_ROOT3_DBL_MAX) {
  127.     OVERFLOW_ERROR(result);
  128.   }
  129.   else {
  130.     result->val = M_PI/(2.0*x) * (1.0 + 3.0/x * (1.0 + 1.0/x));
  131.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  132.     CHECK_UNDERFLOW(result);
  133.     return GSL_SUCCESS;
  134.   }
  135. }
  136.  
  137.  
  138. int gsl_sf_bessel_kl_scaled_e(int l, const double x, gsl_sf_result * result)
  139. {
  140.   if(l < 0 || x <= 0.0) {
  141.     DOMAIN_ERROR(result);
  142.   }
  143.   else if(l == 0) {
  144.     return gsl_sf_bessel_k0_scaled_e(x, result);
  145.   }
  146.   else if(l == 1) {
  147.     return gsl_sf_bessel_k1_scaled_e(x, result);
  148.   }
  149.   else if(l == 2) {
  150.     return gsl_sf_bessel_k2_scaled_e(x, result);
  151.   }
  152.   else if(x < 3.0) {
  153.     return bessel_kl_scaled_small_x(l, x, result);
  154.   }
  155.   else if(GSL_ROOT3_DBL_EPSILON * x > (l*l + l + 1)) {
  156.     int status = gsl_sf_bessel_Knu_scaled_asympx_e(l + 0.5, x, result);
  157.     double pre = sqrt((0.5*M_PI)/x);
  158.     result->val *= pre;
  159.     result->err *= pre;
  160.     return status;
  161.   }
  162.   else if(GSL_MIN(0.29/(l*l+1.0), 0.5/(l*l+1.0+x*x)) < GSL_ROOT3_DBL_EPSILON) {
  163.     int status = gsl_sf_bessel_Knu_scaled_asymp_unif_e(l + 0.5, x, result);
  164.     double pre = sqrt((0.5*M_PI)/x);
  165.     result->val *= pre;
  166.     result->err *= pre;
  167.     return status;
  168.   }
  169.   else {
  170.     /* recurse upward */
  171.     gsl_sf_result r_bk;
  172.     gsl_sf_result r_bkm;
  173.     int stat_1 = gsl_sf_bessel_k1_scaled_e(x, &r_bk);
  174.     int stat_0 = gsl_sf_bessel_k0_scaled_e(x, &r_bkm);
  175.     double bkp;
  176.     double bk  = r_bk.val;
  177.     double bkm = r_bkm.val;
  178.     int j;
  179.     for(j=1; j<l; j++) { 
  180.       bkp = (2*j+1)/x*bk + bkm;
  181.       bkm = bk;
  182.       bk  = bkp;
  183.     }
  184.     result->val  = bk;
  185.     result->err  = fabs(bk) * (fabs(r_bk.err/r_bk.val) + fabs(r_bkm.err/r_bkm.val));
  186.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  187.  
  188.     return GSL_ERROR_SELECT_2(stat_1, stat_0);
  189.   }
  190. }
  191.  
  192. int gsl_sf_bessel_kl_scaled_array(const int lmax, const double x, double * result_array)
  193. {
  194.   if(lmax < 1 || x <= 0.0) {
  195.     GSL_ERROR("domain error", GSL_EDOM);
  196.   }
  197.   else {
  198.     int ell;
  199.     double kellp1, kell, kellm1;
  200.     gsl_sf_result r_kell;
  201.     gsl_sf_result r_kellm1;
  202.     gsl_sf_bessel_k1_scaled_e(x, &r_kell);
  203.     gsl_sf_bessel_k0_scaled_e(x, &r_kellm1);
  204.     kell   = r_kell.val;
  205.     kellm1 = r_kellm1.val;
  206.     result_array[0] = kellm1;
  207.     result_array[1] = kell;
  208.     for(ell = 1; ell < lmax; ell++) {
  209.       kellp1 = (2*ell+1)/x * kell + kellm1;
  210.       result_array[ell+1] = kellp1;
  211.       kellm1 = kell;
  212.       kell   = kellp1;
  213.     }
  214.     return GSL_SUCCESS;
  215.   }
  216. }
  217.  
  218.  
  219. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  220.  
  221. #include "eval.h"
  222.  
  223. double gsl_sf_bessel_k0_scaled(const double x)
  224. {
  225.   EVAL_RESULT(gsl_sf_bessel_k0_scaled_e(x, &result));
  226. }
  227.  
  228. double gsl_sf_bessel_k1_scaled(const double x)
  229. {
  230.   EVAL_RESULT(gsl_sf_bessel_k1_scaled_e(x, &result));
  231. }
  232.  
  233. double gsl_sf_bessel_k2_scaled(const double x)
  234. {
  235.   EVAL_RESULT(gsl_sf_bessel_k2_scaled_e(x, &result));
  236. }
  237.  
  238. double gsl_sf_bessel_kl_scaled(const int l, const double x)
  239. {
  240.   EVAL_RESULT(gsl_sf_bessel_kl_scaled_e(l, x, &result));
  241. }
  242.  
  243.  
  244.